The Jamming Transition in Granular Systems 
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Recent simulations have predicted that near jamming for collections of spherical particles, there will be a 
discontinuous increase in the mean contact number, Z, at a critical volume fraction, <|) c , Above § c , Z and the 
pressure, P, are predicted to increase as power laws in (|> — <|> c . In experiments using photoelastic disks we 
corroborate a rapid increase in Z at (|) e and power-law behavior above (|) c for Z and P. Specifically we find 
power-law increase as a function of — <|) c for Z — Z c with an exponent (3 around 0.5, and for P with an exponent 
\|/ around 1.1. These exponents are in good agreement with simulations. We also find reasonable agreement 
with a recent mean-field theory for frictionless particles. 

PACS numbers: 64.60.-i,83.80.Fg,45.70.-n 
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A solid, in contrast to a fluid, is characterized by me- 
chanical stability that implies a finite resistance to shear and 
isotropic deformation. While such stability can originate from 
long-range crystalline order, there is no general agreement on 
how mechanical stability arises for disordered systems, such 
as molecular and colloidal glasses, gels, foams, and granu- 
lar packings [1]. For a granular system in particular, a key 
question concerns how stability occurs when the packing frac- 
tion, (j), increases from below to above a critical value (j) 6 . for 
which there are just enough contacts per particle, Z, to sat- 
isfy the conditions of mechanical stability. In recent simula- 
tions on frictionless systems it was found that Z exhibits a dis- 
continuity at ()) c followed by a power law increase for (j) > (j) t 
ilia. The pressure is also predicted to increase as a 
power-law above (|) c . 

A number of recent theoretical studies address jamming, 
and we note work that may be relevant to granular systems. 
Silbert, O'Hern et al. have shown in computer simulations 
of frictionless particles J2, S 0] that: a) for increasing (j), 
Z increases discontinuously at the transition point from zero 
to a finite number, Z c , corresponding to the isostatic value 
(needed for mechanical stability); b) for both two- and three- 
dimensional systems, Z is expected to continue increasing as 
((]) — (j) c )P above (j) c , where (5 = 0.5; c) the pressure, P, is ex- 
pected to grow above (|) c as ((j) — § c )^, where \|/ = a/ — 1 in 
the simulations, and a/ is the exponent for the interparticle 
potential. More recent simulations by Donev et. al. for hard 
spheres in three dimensions found a slightly higher value for 
p, p ~ 0.6, in maximally random jammed packings ||5J. It is 
interesting to note that a model for foam exhibits quite similar 
behavior for Z |fj. Henkes and Chakraborty [7] constructed 
a mean field theory of the jamming transition in 2D based on 
entropy arguments. These authors predict power-law scaling 
for P and Z in terms of a variable a, which is the pressure 
derivative of the entropy. By eliminating a, one obtains an al- 
gebraic relation between P and Z — Z t from these predictions, 
which we present below in the context of our data. 

While the simulations agree among themselves at least 
qualitatively, so far, these novel features have not been iden- 



tified in experiments. Hence, it is crucial to test these pre- 
dictions experimentally. In the following, we present experi- 
mental data for Z and P vs. (j), based on a method that yields 
accurate determination of the of contacts and identifies power 
laws in Z and P for a two-dimensional experimental system of 
photoelastic disks. By measuring both P and Z, we can also 
obtain a sharper value for the critical packing fraction (j) c , for 
the onset of jamming, and we can test the model of Henkes 
and Chakraborty. 

The relevant simulations have been carried out predomi- 
nantly for frictionless particles. For real frictional particles 
there will clearly be some differences. For instance, in the 
isostatic limit, Z equals 4 for frictionless disks, whereas for 
frictional disks, Z is around 3, depending on the system de- 
tails 1 8]. Other predictions such as specific critical exponents 
may also need modification. However, one might hope that 
the observed experimental behavior, in particular critical ex- 
ponents, might be similar to that for frictionless particles if 
the frictional forces are typically small relative to the normal 
forces. Indeed, in recent experiments, the typical inter-grain 
frictional forces in a physical granular system were found to 
be only about 10% of the normal forces (9J]. 

Fig. 0J shows a schematic of the apparatus. We use a 
bidisperse mixture (80% small and 20% large particles) of 
approximately 3000 polymer (PSM-4) photoelastic (birefrin- 
gent under stress) disks with diameter 0.74 cm or 0.86 cm. 
This ratio preserves a disordered system. The disks have 
Young's modulus of 4 MPa, and a static coefficient of friction 
of 0.85. The model granular system is confined in a biaxial 
test cell (42cm x 42cm with two movable walls) which rests 
on a smooth Plexiglas sheet. The displacements of the walls 
can be set very precisely with stepper motors. The linear dis- 
placement step size used in this experiment is 40/jm, which 
is approximately 0.005Z), where D is the average diameter of 
the disks. The deformation 8 per particle is less than 1 % in the 
compressed state. The setup is horizontal and placed between 
crossed circular polarizers. It is imaged from above with an 8 
MP CCD color camera which captures roughly 1200 disks in 
the center of the cell, enabling us to visualize the stress field 
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FIG. 1 : (a) Schematic cross-section of biaxial cell experiment (not 
to scale). Two walls can be moved independently to obtain a de- 
sired sample deformation, (b) Examples of contacts and particles 
that are either close but not actually in contact, or contacts with very 
small forces. Circles show true contacts, squares show false apparent 
contacts, (c) Image of a single disk at the typical resolution of the 
experiment, (d) Sample image of highly jammed/compressed state 
and (e) almost unjammed state. 



within each disk (Fig.Q]). We then obtain good measurements 
of the vector contact forces (normal and tangential = frictional 
components) 131 . 

We also use the particle photoelasticity to accurately deter- 
mine the presence or absence of contacts between particles. 
In numerical studies one can use a simple overlap criterion 
to determine contacts: a contact occurs if the distance be- 
tween particle centers is smaller than the sum of the particle 
radii. However, in experimental systems, a criterion based 
solely on the particle centers is susceptible to relatively large 
errors which include false positives (Fig.QJ) - squares) as well 
as false negatives (circles). As seen in Fig. QJ), the contacts 
through which there is force transmission appear as source 
points for the stress pattern. Further details are given in the 
supplementary material, section I. 

We use two protocols to produce different packing frac- 
tions: we either compress the system from an initially stress- 
free state, or decompress the system until the end state is es- 
sentially a stress-free state. The results for both protocols are 
the same within error bars above (j) c ; below jamming, the data 
for Z obtained by compression are a few percent below those 
for decompression. Below, we will present decompression 
data. Figures [T}l,e show the initial highly stressed state and 
the end state after decompression, respectively. After each 
decompression step, we apply tapping to relax stress in the 
system. This could be seen as roughly analogous to the an- 
nealing process invoked in some simulations. Two images are 
captured at each state: one without polarizers to determine the 
disk centers and one with polarizers to record the stress. 
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FIG. 2: Average contact number and pressure at the jamming tran- 
sition. Top and bottom panels show Z — Z c and P vs <|> — <|) f , respec- 
tively, with rattlers included (stars) or excluded (diamonds). Dashed 
and full curves in the top panel give power-law fits (0 — <|) C )P with 
P = 0.495 and 0.561 for the case with and without rattlers, respec- 
tively. Full curve in the lower panel gives the fit (0 — <^ c )^ with 
\j/ = 1.1; dashed line shows a linear law for comparison. Inset: Z 
vs for a larger range in 0. 



The average Z can be computed either by counting only the 
force bearing disks or by counting all the disks including rat- 
tlers which do not contribute to the mechanical stability of the 
system. We consider as rattlers, all the disks which have less 
than 2 contacts. For the number of rattlers beyond the transi- 
tion point we find an exponential decrease with — C ; hence, 
a divergence in the number of rattlers at (j) c is not indicated by 
the data. 

We next compute the Cauchy stress tensor for each disk, 
Ojj = jfiY<(FiXj + FjXj); P, is the trace of this tensor. Here, 
A is the Voronoi area for the given disk, and the sum is taken 
over contacts for a given disk. We then compute the average of 
the pressure over the ensemble of disks in the system. For the 
data presented below, we performed two sets of experiments: 
one with a larger range, 0.8390 < < 0.8650, and also larger 
step size, A(j) = 0.016, and - after the jamming region was 
identified - a second set at a finer scale with 0.840745 < (j) < 
0.853312, with a step size, A(j) = 0.000324. 
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The inset in Fig. [2] shows data for Z over a broad range of 
(j) (with rattlers-stars; without-squares). These data show a 
significant rise in Z at the jamming transition. While this rise 
is not sharply discontinuous, it occurs over a very small range 
in (j). At higher (]), the variations of the curves are similar with 
and without rattlers. At lower (j), their behavior differs: The 
values of Z drop lower for the case with rattlers. The pressure 
P((|) — (j) 6 ) in Fig. [2] shows a fiat background below jamming, 
and then a sharp positive change in slope at a well defined (]). 
The pressure is not identically zero below jamming for similar 
reasons that the jump in Z is not perfectly sharp, as discussed 
below. 

To compare these experimental results to predictions above 
(|) c , we carry out least squares fits of Z — Z c and P to (j) — (j) c . 
These fits depend on the choice of § c , which has some ambi- 
guity due to the rounding; the data allow a range from around 
0.840 to 0.843. In fact, (|) c can be determined in several ways: 
the point where Z reaches 3, the point where P begins to rise 
above the background, etc. (cf. supplemetary material). We 
show results of these fits in Fig. [2] starting with the upper 
panel, which shows power-law fits (Z — Z c ) °^ ((j) — (jv)' 3 . The 
fitted exponent p depends on the choice of <]) e but the variation 
is small without rattlers, 0.494 < p < 0.564, and somewhat 
larger with rattlers, 0.363 < p < 0.525. The details for sev- 
eral different specific fits are given in the supplementary ma- 
terial, section II. The point (|) c = 0.84220 where P rises above 
the background level is used in Fig. [2] and yields a consis- 
tent fit for both P and Z. The point where Z reaches 3 for the 
case without rattlers agrees with the previous case to within 
8(j) c - = 0.0005, and the exponents are quite similar. Comparing 
with the simulations for frictionless particles, we find that our 
values of P « 0.55 for the data without rattlers are larger than 
the value of 0.5 reported in , but smaller than those of 
Donev et. al. |5|] obtaining 0.6 in 3D. In contrast, for a model 
of frictional disks under shear, Aharonov and Sparks ifioll ob- 
tain the much lower value of 0.36. However, a direct com- 
parison is not possible to the present case of jamming under 
isotropic conditions. 

Figure |2] shows the variation of P with (]) in the lower panel, 
indicating a clear transition at C — 0.8422 ±0.0005. For this 
choice of <|) c , P increases as foe (<|)_(]) c ) v l' with \\t — 1.1 ±0.05 
above (|) c . This value of \|/ pertains to a fit over the full range 
(j) > (j) c of Fig. |2j a larger exponent would be obtained if the fit 
range were limited to very close to § c . This value is close to 
the value \\i = 1 .0 found ||2|, |3|] for a linear force law, and this 
linear law is indicated as a dashed line in Fig. [2] One expects 
such a linear force law (with a logarithmic correction) for ideal 
disks, but direct mechanical calibration of the force law for the 
cylinders is closer to 8 3 / 2 (see supplementary material). This 
rather high exponent for the force law is attributable to the 
small asperities, which influence the force law for small de- 
formations. However, the photoelastic response is detectable 
only for 8 > 150/jm, and for such 8's, the force law is close to 
locally linear in 8. 

From the P vs. (j) data, we can also obtain the bulk modulus, 
B = —AdP/dA, where A is the area enclosed by the system 
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FIG. 3: Results from new computer simulations. For all fits § c = 
0.84005. (a) Average overlap per particle in units of the mean particle 
radius is linear in <)) — <|) c . (b) P obtained from the Cauchy stress 
tensor (circles) and the force on the walls (squares) satisfy a power 
law ((]) — <))c) with \\l = 1.13; dashed line shows a linear law for 
comparison, (c) Z (rattlers included) exhibits a power law Z — Z c oc 
(<|) - tyc f with Z c = 3.94 and (3 = 0.5015. 



boundaries. Since, (j) = A p /A, where A p is the (presumably 
fixed) area occupied by the disks, B = §dP/d§. Then, B ^ 
(ty-tyc)^- 1 , which gives a weak pressure variation of B above 
(j) c . We note that anomalous results for the bulk modulus have 
been observed in acoustical experiments by Jia, and discussed 
by Makse et. al. [11], where the bulk modulus near (j) 6 varied 
faster with P than was previously expected because of changes 
in Z. 

Since P in Fig. [2] corresponds closely to expectations for a 
linear force law, we performed a computer simulation for a 
poly disperse system of 1950 particles with a linear force law 
(k n = 10 5 N/m) without friction; details can be found in 11211 . 
In Fig. [3] the results are shown for a larger range in density 
than done in earlier studies. All the data in Fig.[3]can be fitted 
with a single value for the transition density of (j) 6 = 0.84005. 
While the average overlap per particle (equivalent of the de- 
formation 8 for physical particles) is clearly linear in (]), the 
pressure P is not: P increases faster than linear with an ex- 
ponent close to the one found in the experiment. Z is also 
consistent with a power-law exponent close to 0.5. With the 
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FIG. 4: Pressure vs. Z; Experimental data and a fit to the model of 
Henkes and Chakraborty [7]. In this fit, the constant C defined in the 
text is treated as an adjustable parameter. The other fitting parameter 
isZ c . 

rattlers included, Z at (|) c , Z c — 3.94, is slightly below the iso- 
static value of 4 for a frictionless system of disks. 

To connect with the predictions of Henkes and Chakraborty 
J3l, we consider P — P c vs. Z. The prediction from their 
Eq. (10) is equivalent to (P - P c ) jP c =u- [(4m 2 + l) 1 / 2 - 
l]/2, where u = C(Z — Z c ) and C — e/a c - is a system- 
dependent constant. Thus, e is a measure of the grain elas- 
ticity, and e = corresponds to completely rigid grains. Also, 
a c - is the critical value for a. In fitting to this form, we may 
adjust P c , (within reason) Z c , and C. In Fig. |4] we find rea- 
sonable although not perfect agreement with this prediction 
(above (j) c ), and obtain Z c = 3.04, which is close to the iso- 
static value Z c = 3. 

We now turn to the rounding that we observe in Z quite 
close to the transition, and the background pressure that we 
obtain near (j) c . One possible explanation is the friction be- 
tween the disks and the Plexiglas base. This could help freeze 
in contact forces and contacts. However, a simple estimate 
of the upper bound for the friction with the base shows that 
this cannot be a significant effect, at least as regards the pres- 
sure background. To obtain an estimated upper bound for 
the base friction on P, we assume that base friction can sup- 
port inter-grain contact forces corresponding to the maximum 
base frictional force per grain, Ff = /Jba>ng = 2.8 x 10~ 3 N, 
where m is the mass of a grain and \i\, a < 1 is the friction be- 
tween a particle and the base. Assuming Z inter-particle con- 
tacts and one particle-base contact per grain, we estimate the 
resulting upper bound on the perturbation to the pressure as 
5P ~ (ZF f R)/{nR 2 ) ~ 0.22Z N/m, where R is a disk radius. 
Since Z ~ 3, this pressure is almost two orders of magnitude 
too small to be of relevance. An additional issue concerns the 
anisotropy that is induced during compression or expansion 
by the apparatus. This induced anisotropy is difficult to avoid 
and/or relax close to (|) c even in the simulation, but it remains 



small. It is visible in Fig. QJ, where a weak array of force 
chains tends to slant from lower left to upper right. Among 
other reasons, the anisotropy can be induced by wall friction 
due to the confining lateral boundaries of the biaxial appara- 
tus. 

We conclude by noting that these experiments, the first of 
which we are aware, demonstrate the critical nature of jam- 
ming in a real granular material. Our results take advantage of 
the high accuracy in contact number Z that is afforded when 
the particles are photoelastic. Z shows a very rapid rise at a 
packing density (j) c = 0.8422. The fine resolution in density 
allows us to see that the transition is not as sharply discontin- 
uous under the present experimental conditions as in the com- 
puter simulation. Above <|) c , Z and P follow power laws in 
(]) — <|) c with respective exponents p of 0.5 to 0.6 and \\t ks 1.1. 
The values for both p and \\i are consistent with recent sim- 
ulation results for frictionless particles. In addition, we find 
reasonable agreement with a mean field model of the granu- 
lar jamming transition, again for frictionless particles. These 
results suggest that effects of friction on jamming are likely 
modest, although perhaps not ignorable. That jamming in the 
experiment occurs over a narrow, but finite range in (j) seems 
mostly to be caused by small residual shear stresses that are 
induced by interactions with the walls confining the sample 
(not the base supporting the particles). The ability of a small 
amount of shear to affect the jamming transition is interest- 
ing, and points to the need for a deeper understanding of the 
effects of anisotropy. 
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I. SUPPLEMENTARY METHODS 

In this supplement, we provide experimental details which 
we discuss in the context of Fig. 1. A key point concerning 
the experiments is the use of photoelasticity (stress-induced 
birefringence) to obtain vector forces at interparticle contacts. 
This technique has the added advantage of determining to 
good accuracy whether a contact is present or not. 

a. Photoelastic Method and Determination of Contacts 
A stressed photoelastic particle (in our case, a disk) when 
viewed through crossed circular polarizers, shows a pattern 
of light and dark bands. The light rays traversing the polariz- 
ers and a particle (along the axial direction of the disk) have 
an intensity / = I sin 2 [{<3\ — 02)C\. Here, the a, are the prin- 
ciple stresses within the particle; C is a constant that depends 
on the the thickness and properties of the disk, and on the 
wavelength of the light 01- Given a set of contacts for a disk, 
and forces at these contacts, the specific photoelastic pattern 
is determined. Here, we take advantage of the fact that a two- 
dimensional description for the stresses is appropriate. As- 
suming that the contact forces are well-approximated as point- 
like, the Boussinesq solution gives the stresses within the disk 
For these experiments, we solve the inverse problem: we 
have the light intensities of the photoelastic pattern within a 
disk, and we find the contact forces. We use an automated 
computer algorithm which uses the vector contact forces as 
nonlinear least-squares fit parameters. The fitting procedure 
minimizes differences between the experimentally measured 
intensity pattern for a disk and the intensity pattern that would 
be obtained for a given set of contact forces 

In order to improve the discrimination between false and 
true contacts we employ a two step process. The first step in- 
volves obtaining possible contacts based on the distances be- 
tween disk centers; if the particle centers are within D ± 0. ID, 
where D is the mean center-to-center distance of a particle 
pair, the disks are considered to be in potential contact. This 
estimate of contacts is markedly improved by utilizing the 
photoelastic stress images at various exposure times for each 
state, such that eventually most of the force transmitting con- 
tacts can be seen. As seen in Fig. lb, the contacts through 
which there is force transmission, appear as source points 
for the stress pattern. This effect can be quantified by mea- 
suring the intensity and the gradient square of the intensity 
(G 2 = |V/| 2 where the gradient is taken in the plane of the 
disk) around the contact 14]]. A true, force bearing contact 
can be distinguished by employing appropriate thresholds in 
intensity and in G 2 . The thresholds in intensity and in G 2 are 
useful in capturing contacts with very small forces, since these 
quantities are higher near force bearing contacts. The final er- 



ror in average Z is around 3.5% for rather low (]), and around 
1.5% for higher (j). 

b. Calibration of the Force Law A direct mechanical 
calibration for the particles using a digital force gage is shown 
in Fig. IS II The dotted curve shows a force law F °^ 8 3 / 2 . A 
linear fit describes the calibration data well for 8 > 250/jm 
which is comparable to the surface roughness of the cylinders. 
The photoelastic response is detectable for displacements that 
exceed the right end of the gray bar at 8 rj 150/jm. In the ef- 
fective range for the photoelastic technique, 8 > 150/jm, the 
force vs. displacement curve is reasonably well described by 
a straight line. 




0.3 
8(mm) 

FIG. SI: Calibration of the contact force F for a representative disk 
pushed against a hard surface by a displacement 8. The experimental 
data (squares) are fitted by the power law F = 2.52N (S 1,54 ) (dotted) 
and by the linear law F = 2.56N (8 — 0.16) (full curve). Here, all 
lengths are given in mm. The gray bar indicates the roughness of the 
cylinder surface. Photoelastic response is reliably detectable to the 
right of this bar. 



H. SUPPLEMENTARY TABLE 

c. Details of the Fitting Procedure For the fit of the data 
with the power law Z — Z c = a(0 — C )P we examine a range 
of values for <|) c and obtain the exponents for the power-law 
fits given in Table [Si] Here, (j) c is selected, and Z c , and p are 
the fitting parameters. For the case without rattlers, p ranges 
from 0.49 to 0.56, and Z e ranges from 2.40 to 3.08. For the 
case with rattlers, p shows more variation (0.36 - 0.52), and 
the errors in Z c are larger. For the entire range of (j), the root 
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<t>c 


Without Rattlers 


With Rattlers 




z c 


P 


RMSE 


Z c 


P 


RMSE 


0.84058 


2.397±0.135 


0.517±0.064 


0.049 


1.198±0.310 


0.502 ±0.093 


0.109 


0.84075 


2.512±0.138 


0.547 ±0.073 


0.051 


1.071 ±0.359 


0.460 ±0.090 


0.103 


0.84172 


2.632±0.151 


0.494 ±0.077 


0.045 


0.9747 ±0.458 


0.363 ±0.083 


0.080 


0.84204 


2.858 ±0.127 


0.564 ±0.086 


0.045 


1.183 ±0.413 


0.367 ±0.079 


0.072 


0.84220 


2.838±0.171 


0.533 ±0.102 


0.046 


1.490 ±0.427 


0.405 ±0.096 


0.072 


0.84236 


2.916±0.133 


0.556 ±0.093 


0.046 


1.744 ±0.298 


0.445 ±0.088 


0.075 


0.84269 


3.003 ±0.124 


0.563 ±0.095 


0.043 


1.989±0.267 


0.469 ±0.092 


0.071 


0.84301 


3.075 ±0.120 


0.560 ±0.095 


0.041 


2.280 ±0.235 


0.525 ±0.108 


0.072 



TABLE SI: Power-law exponents and critical contact numbers obtained as fitting parameters, at various critical packing fractions. The RMSE 
gives the root mean squared errors for the fits. The indicated uncertainties in both Z c , and (3 are obtained from the 95% confidence interval of 
the best-fit parameter values. 



mean squared errors (RMSE) are larger for the case with rat- - 0.051), indicating that power-law fits are consistently better 
tiers (0.071 - 0.109), than for the case without rattlers (0.041 when rattlers are excluded. 
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